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STOCHASTIC MODELING IN NANOSCALE BIOPHYSICS: 
SUBDIFFUSION WITHIN PROTEINS 

By S. C. Kou 1 

Harvard University 

Advances in nanotechnology have allowed scientists to study bio- 
logical processes on an unprecedented nanoscale molecule-by-molecule 
basis, opening the door to addressing many important biological 
problems. A phenomenon observed in recent nanoscale single-molecule 
biophysics experiments is subdiffusion, which largely departs from 
the classical Brownian diffusion theory. In this paper, by incorporat- 
ing fractional Gaussian noise into the generalized Langevin equation, 
we formulate a model to describe subdiffusion. We conduct a de- 
tailed analysis of the model, including (i) a spectral analysis of the 
stochastic integro-differential equations introduced in the model and 
(ii) a microscopic derivation of the model from a system of interacting 
particles. In addition to its analytical tractability and clear physical 
underpinning, the model is capable of explaining data collected in 
fluorescence studies on single protein molecules. Excellent agreement 
between the model prediction and the single-molecule experimental 
data is seen. 

1. Introduction. 

1.1. Background: Nanoscale single- molecule biophysics. It is said that 
the famous Richard Feynman once described seeing the images of single 
atoms as a "religious experience" for physicists. Recent advances in nan- 
otechnology have turned Feynman's "religious" encounter into daily reality. 
In particular, scientists are now able to study biological processes on an 
unprecedented nanoscale molecule-by-molecule basis [cf. Moerner (2002), 
Nie and Zare (1997), Tamarat, Maali, Lounis and Orrit (2000), Weiss (2000), 
Xie and Lu (1999), Xie and Trautman (1998), Kou, Xie and Liu (2005)], thus 
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opening the door to addressing many problems that were inaccessible just a 
few decades ago. 

Compared with traditional experiments, which involve a population of 
molecules, (nanoscale) single-molecule experiments offer many advantages. 
First, they provide experimental data with more accuracy and higher res- 
olution because scientists can "zoom in" on individual molecules to study 
and measure them. Second, by following individual molecules, these single- 
molecule experiments can capture transient intermediates and detailed dy- 
namics of biological processes. This type of information is rarely available 
from the traditional population experiments. Third, in a living cell, many 
important biological functions are often carried out by single molecules; thus, 
understanding the behavior of molecules at the individual level is of crucial 
importance. Many new discoveries [see, e.g., Asbury, Fehr and Block (2003), 
English et al. (2006), Kou et al. (2005), Lu, Xun and Xie (1998), Zhuang et 
al. (2002)] have emerged from the nanoscale single-molecule studies. 

The technological advance also brings opportunities and challenges for 
stochastic modelers because the individual molecules, subject to statisti- 
cal and quantum mechanics of the nanometer world, behave stochastically. 
Characterizing their fluctuation thus requires stochastic models [Kou (2007)]. 
In the current article we focus on modeling the phenomenon of subdiffusion 
observed in single-molecule experiments to exemplify the stochastic model- 
ing problems in the field. 

1.2. Subdiffusion in proteins: The experimental finding. Since Einstein's 
and Wiener's ground breaking works in the early 20th century, the the- 
ory of Brownian motion and diffusion processes has revolutionized not only 
physics, chemistry and biology, but also probability and statistics. One 
key characteristic of Brownian motion is that the second moment E[x 2 (t)], 
which in physics corresponds to the mean squared displacement (location) 
of a Brownian particle, is proportional to time t. In some systems [cf. 
Bouchaud and Georges (1990), Klafter, Shlesinger and Zumofen (1996), 
Sokolov, Klafter and Blumen (2002)], scientists, however, have discovered 
a clear departure from Brownian diffusion. The mean squared displacement 
E[x 2 (t)] there is no longer proportional to t, but rather E[x 2 (t)] cx t a , where 
< a < 1. Because a < 1, these movements satisfying E[x 2 (t)] oc t a are de- 
fined as subdiffusion. Recent single-molecule biophysics experiments [Yang 
et al. (2003), Kou and Xie (2004), Min et al. (2005)] reveal that subdiffusion 
may be quite prevalent in biological systems. 

In a 2003 Science paper [Yang et al. (2003)] scientists conducting single- 
molecule experiments on a protein-enzyme system observed this subdiffusion 
phenomenon. The experiment studied a protein-enzyme compound, called 
Fre, which is involved in the DNA synthesis of the bacterium E. Coli. In 
the reactions Fre works as a catalyst. Figure 1 shows the crystal structure 
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Fig. 1. The crystal structure of Fre. The two substructures FAD and Tyr are shown. 

of Fre, which contains two smaller structures: FAD (an electron carrier) and 
Tyr (an amino acid). The 3D conformation (shape) of Fre spontaneously 
fluctuates, and consequently, the (edge-to-edge) distance between the two 
substructures FAD and Tyr varies over time. It was found in the experiment 
that the stochastic distance fluctuation between FAD and Tyr undergoes a 
subdiffusion. Section 5 provides more details about the experiment. 

1.3. Modeling subdiffusion. To explain this subdiffusion phenomenon, 
we shall formulate a stochastic model by incorporating fractional Gaus- 
sian noise (formally the derivative of fractional Brownian motion) into a 
stochastic integro-differential equation framework governed by the general- 
ized Langevin equation. 

Since its introduction [Mandelbrot and Van Ness (1968)], fractional Brow- 
nian motion (fBm) has proven to be an indispensable tool for stochastic 
modeling [Samorodnitsky and Taqqu (1994), Taqqu (1986), Whitt (2002), 
Adler, Feldman and Taqqu (1998)]; its applications range from 
queuing systems [Glynn and Zeevi (2000), Heath, Resnick and Samorodnitsky 
(1997), Konstantopoulos and Lin (1996)] to finance [Heyde (1999), Mandelbrot 
(1997)] to internet traffic [Leland, Taqqu, Willinger and Wilson (1994), 
Crovella and Bestavros (1996), Mikosch, Resnick, Rootzen and Stegeman 
(2002)]. On the other hand, the generalized Langevin equation [Chandler 
(1987), Zwanzig (2001), Wang and Tokuyama (1999)], used primarily in 
the physics literature, has attracted less attention from probabilists and 
statisticians. Notably, it is the connection between fBm and the generalized 
Langevin equation (GLE), as we shall see, that leads to a satisfactory model 
to account for the experimentally observed subdiffusion within proteins. 
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A key requirement in the construction of biophysical models, in addition 
to the preference of analytical tractability, is that the model must agree with 
fundamental physical laws and should have a sound physical foundation. 
Hence, we will also discuss the model's physical basis. 

From an applied probabilistic/statistical point of view, to describe sub- 
diffusion, we study several stochastic integro-differential equations driven 
by fBm. We conduct a detailed analysis, in particular, a spectral analysis, 
of their properties (in Sections 2 and 3). We also consider the connection 
between our model and interacting particle systems through a microscopic 
derivation of the model from a system of interacting particles (in Section 4) . 
Some of the mathematical structures of GLE and fBm have been indepen- 
dently considered in Kupferman (2004). In this paper, in addition to the 
detailed stochastic investigation, we apply the analytical results to fit the 
nanoscale single-molecule experimental data (in Section 5) and to show that 
the model successfully explains the experimental observations. 

The paper is organized as follows. Section 2 concerns the basic model 
to describe the subdiffusive movement of a free particle. Section 3 studies 
subdiffusive motion under the presence of an outside potential. Section 4 in- 
vestigates the physical foundation of our model. Section 5 applies the model 
to explain the nanoscale single-molecule experimental results, which shows 
close agreement between the model and the data. Section 6 concludes the 
paper with a discussion and raises some open problems. 

2. Modeling subdiffusion of a free particle. 

2.1. Physical Brownian motion via the Langevin equation with white noise. 
To facilitate the discussion of our model, let us first review how the law of 
Brownian diffusion was derived in physics because the Brownian motion 
used by physicists and the term Brownian motion used in probability and 
statistics refer to different things: Physicists' Brownian motion corresponds 
to the integral of the Ornstein-Uhlenbeck process, as we shall see shortly, 
whereas statisticians and probabilists' Brownian motion refers to the Wiener 
process, although both share the characteristic of E[x 2 (t)] oct for large t. 

Suppose we have a Brownian particle with mass m suspended in wa- 
ter. The physical description of the particle's free motion starts from the 
Langevin equation [Risken (1989), Van Kampen (2001), Karlin and Taylor 



where v(t) is the velocity of the particle at time t, and dv{t)/dt is the 
acceleration of the particle. On the right-hand side, £ is the friction constant, 
reflecting the fact that the resistance the particle receives is proportional to 



(1981)] 




m 



dv(t) 
dt 



tv(t) + F(t), 
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its velocity, and F(t) is the white noise, formally the derivative of the Wiener 
process. 

The Langevin equation has an important physical constraint that links 
the friction constant £ with the noise level, because both the movement of 
the particle and the friction originate from one source — the collision between 
the particle and surrounding water molecules. Borrowing physicists' notation 
5(-) of Dirac's delta function, this link can be expressed as 

(2.2) E[F(t)F(s)}=2(k B T-5(t-s), 

where k B is the Boltzmann constant, and T is the underlying temperature. 
This proportional relationship between the noise level and the friction con- 
stant is a consequence of the fluctuation-dissipation theorem in statistical 
mechanics [Chandler (1987), Hill (1986)]. In the more familiar probability 
notation, equations (2.1) and (2.2) translate to 



(2.3) mdv(t) = -Cv(t)dt + y/2Ck B TdB(t), 

where B(t) is the Wiener process, and the formal association of u F(t) = 
y/2(k B TdB(t)/dt" is recognized. 

The stationary solution of equation (2.3) is the Ornstein-Uhlenbeck pro- 
cess [Karlin and Taylor (1981)], which is Gaussian with mean function 
E[v(t)] = and covariance function E[v(t)v(s)] = ^^exp(— £-\t — s\). It 
follows that, for the displacement, x(t) = Jq u(s) ds, which is the actual ob- 
served motion, the variance is 

Var[a?(t)] =E[x 2 {t)} = f f E[v{s)v{u)]duds 
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c 2 



~ 2 B t, for large t. 

The last line is known (in physics) as Einstein's Brownian diffusion law 
(where 2k B T/Q is the Einstein diffusion constant). It is worth emphasizing 
that when physicists talk about the Brownian motion of a free particle, they 
refer to the integral of the Ornstein-Uhlenbeck process and the correspond- 
ing equation (2.4), which bears the resemblance of Var[a;(t)] cx t for large t 
to the Wiener process. 



2.2. Toward subdiffusion. The classical theory of Brownian diffusion, 
however, cannot explain the subdiffusion phenomenon, which, defined by 
Var[x(t)] oc t a with < a < 1 for large t, has been notably observed in dis- 
tance fluctuation within proteins. We will explain the experimental details 
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in Section 5, where the theoretical results and predictions of our model are 
compared with the experimental data. In this subsection and the next we 
focus on the model itself. 

The starting point of our model is the generalized Langevin equation 
(GLE) [Chandler (1987), Zwanzig (2001)] 

(2.5) m ( ^L = -Q[ t v(u)K{t-u)du + G{t), 

at J-oo 

where, in comparison with the Langevin equation (2.1), (i) a noise G(t) 
having memory replaces the memoryless white noise and (ii) a kernel K 
convoluted with the velocity makes the process non-Markovian. 

The reason that both K and G{t) appear in the equation is that any 
closed (equilibrium) physical system must satisfy the fluctuation-dissipation 
theorem, which requires the memory kernel K(t) and the fluctuating noise 
to be linked by 

(2.6) E[G(t)G(s)]=k B T(-K{t-s). 

In an intuitive sense this relationship arises because both the friction and 
the motion of the particle originate from the collision between the particle 
and its surrounding media. Equations (2.5) and (2.6) contain their Langevin 
counterpart (2.1) and (2.2) as a special case, in which the kernel K is the 
delta function. 

It is important to note that relationship (2.6) also rules out models like 
mdv(t)/dt = —(v(t) + G(t), with G(t) an arbitrary noise to describe biophys- 
ical processes because this kind of model violates the fluctuation-dissipation 
theorem, and thus cannot correspond to any (equilibrium) physical process. 
In other words, to have a physically meaningful description of biophysi- 
cal processes, the subdiffusion process, in particular, the convolution term 
f-oo v ( u )K(t — u) du must be present in equation (2.5). 

Under the framework of GLE, the key question is as follows: Is there a 
combination of kernel function and noise structure that can lead to sub- 
diffusion? To answer this question, we note that the white noise is math- 
ematically interpreted as the formal derivative of the Wiener process. It 
is well known that the Wiener process is the unique process characterized 
by (i) being Gaussian, (ii) having independent increments, (iii) having sta- 
tionary increments, and (iv) being self-similar. These properties carry their 
physical meanings: the independent increments of the Wiener process make 
the white noise independent across time; the stationary increments of the 
Wiener process mean that the white noise is time translation invariant; the 
self-similarity means that the white noise is invariant to time scale change. 
To generalize the white noise, we want to maintain as many nice properties 
as possible and at the same time introduce memory. This leads us to consider 
processes with the following three properties: (i) Gaussian, (ii) stationary 
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increments, and (iii) self-similar. The only class of processes that embod- 
ies all three properties is the fBm process, Bn(t) [Embrechts and Maejima 
(2002), Samorodnitsky and Taqqu (1994)], which is Gaussian with mean 
E[B H it)} = and covariance function E[B H {t)B H (s)} = \{\t\ 2H + \s\ 2H - 
\t — s\ 2H ). H, between and 1, is the Hurst parameter; Bjj{t) reduces to 
the Wiener process when H = 1/2. 

2.3. Subdiffusion via the generalized Langevin equation with fractional 
Gaussian noise. Taking G(t) in (2.5) to be the (formal) derivative of fBm, 
which is referred to as the fractional Gaussian noise (fGn), F}j(t) = 
\Z2^kBTdBjj(t) /dt, we have the following model: 

The model for subdiffusion: 

(2.7) m dt = ~^£ v{u)K H (t-u)du + F H (t), 

where the kernel Knit), according to equation (2.6), is (formally) given by 

K H (t) = E[F H (0)F H (t)]/(k B TC) 

-B H (h) B H (t + h)-B H (ty 



(2i 



21im£ 

HO 



h 



= 2]im^(\t + h\ 2H + \t-h\ 2H -2\t\ 2H ) 
HO 2h 2 

= 2H(2H- l)\t\ 2H ~ 2 , fort t^O. 

In the more familiar probability notation equation (2.7) can be written as 
follows: 

The model for subdiffusion: 

(2.9) mdv(t) = -((J v(u)K H (t - u) duj dt + y / 2(k B T dB H (t), 
where, as in (2.8), 

(2.10) K H {t) = 2H(2H- l)\t\ 2H - 2 , fort t^O. 

One question about this model equation (2.9) arises immediately from 
a probabilistic standpoint: what is the interpretation of the integral with 
respect to fBm? 

Stochastic integrals driven by fBm have been an active research area in 
recent years. The constructions of the stochastic integral include restricting 
the integrand to specific classes of functions as in Dai and Heyde (1996), 
Gripenberg and Norros (1996), Lin (1995), Shiryaev (1998), using path- 
wise integration for the case of 1/2 < H < 1 as in Mikosch and Norvaisa 
(2000), applying Malliavin calculus as in Albs, Mazet and Nualart (2000), 
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Duncan, Hu and Pasik-Duncan (2000), Nualart (2006), and using regular- 
ization as in Rogers (1997), Carmona and Coutin (2000). For reviews, see, 
for example, Duncan, Hu and Pasik-Duncan (2000), Pipiras and Taqqu (2000), 
Pipiras and Taqqu (2001) and Embrechts and Maejima (2002). 

As we shall see shortly, the case ofl/2<if<lis particularly relevant to 
our description of subdiffusion here. In this case, pathwise integration ap- 
pears most natural, and hence, we interpret our model equation (2.9) with 
H > 1/2 in the pathwise Riemann-Stieltjes sense [Mikosch and Norvaisa 
(2000)]. This pathwise integration allows us to treat integrals with respect 
to dB}{{t) as if they were classical integrals, which simplifies our calculation. 

The presence of the convolution term and the dBn{t) term makes equation 
(2.9) non-Markovian and nonstandard. The solution to model (2.9) is given 
by the following theorem, whose proof is deferred to the Appendix. 

Theorem 2.1. Let Kh(u)) and K^{u) denote the Fourier transforms 
of the kernel Knit) on the entire and positive real lines, respectively: 



/oo 
e ltuJ K H (t) dt = 2T(2H + 1) sin(ff7r)|c 
-oo 

K+(u)= e^K H {t)dt 



(2.12) 

= T(2H + l)|w| 1-2H [sin(fl'7r) - zcos(#7r)sign(w)], 

where sign(w) is the sign function, which is 1 ifui>0 and —1 ifu><0. Then 
under the pathwise interpretation of dB}j(t) for 1/2 < H < 1, the solution 
to model (2.9) is 



/oo 
r{t-u) dB H (i 
-oo 



where the deterministic function r(t) is given by 

1 f°° 1 



r(t) = — I e _,<w dw. 

2ir J-oo (Kjj(u!) — imuj 

Furthermore, the solution v(t) is a stationary Gaussian process with mean 
function E[v(t)] =0 and covariance function 

1 r°° 

C v {t) = E[v(0)v(t)} = — / e- ltw C v {u) dw, 
where the Fourier transform C v (ui) of C v (t) is given by 

/oo 
e itul C v {t) dt = k B TCK H (io)/\CK^{uj) - imuj\ 2 . 
-oo 
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Remark 1. When H -> 1/2, we have C v (u) -> 2k B T(/(( 2 + m 2 uj 2 ) 
and S[u(0)u(t)] = C v (t) — ► exp(— -Mi|), which recovers the Ornstein- 
Uhlenbeck result from classical Brownian diffusion. 

Remark 2. A careful reader might note that equations (2.11) and (2.12) 
involve Fourier transforms of power functions, which are not integrable, so a 
natural question here is their meaning. In general, the Fourier transform / of 
a nonintegrable function f is defined as /(w) = lim^^oo Jf^o/(^)exp(-|t|/a) x 
exp(iiw) dt, that is, the limit of exponential damping. It is in this sense that, 
for example, (t)" 1 ^ 2 and -v/^IcjI -1 ^ 2 are regarded as a Fourier transform pair. 
See Champeney (1987) for a thorough discussion. All the Fourier and inverse 
Fourier transforms in this article are defined in this general sense. 

The next theorem, whose proof is given in the Appendix, details how our 
model (2.9) explains subdiffusion. 

Theorem 2.2. Under model (2.9), let x(t) = f^v(s)ds be the displace- 
ment. Then for 1/2 < H < 1, the mean squared displacement 

VMHt)HE[l(t)V M __g^L_ ,- 2H 

oc t 2 ~ 2H , for large t. 
In other words, our model with 1/2 < H < 1 leads to subdiffusion. 

Remark 3. When H — ► 1/2, the right-hand side of (2.14) converges to 
(2ksT ' l X)t, and we recover the Einstein Brownian diffusion law. 

The above theorem says that, for a particle, if its velocity process v(t) is 
governed by model (2.9) with H > 1/2, then the second moment E[x 2 (t)] of 
the displacement x[t) = /q v(s) ds satisfies E[x 2 (t)] oc t 2 ~ 2H for large t, which 
is exactly the characteristic of subdiffusion. Therefore, our model (2.9) with 
Hurst parameter H > 1/2 leads to an explanation of how subdiffusion could 
arise in biophysical systems — as long as the fractional Gaussian noise drives 
the underlying physical fluctuation, the system will be subdiffusive. 

3. Modeling subdiffusion under external potential. The model (2.9) ex- 
plains subdiffusion of a free particle, that is, the motion of a particle without 
the influence of an outside force (or potential). If there is an external po- 
tential U(x) (e.g., a magnetic field), which is a function of the displacement 
x(t), the model has to be modified. 

In the case of white noise and the Langevin equation (2.1), which corre- 
sponds to classical Brownian diffusion, one should add the term —U'(x{t)) 
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to the right-hand side of the equation to account for the external potential 
[Risken (1989), Van Kampen (2001), Karlin and Taylor (1981)]: 

m^P- = -(v(t)-U'(x{t)) + F(t), x(t)= f t v(s)ds. 
at Jo 

For GLE, similarly, the term —U'(x(t)) also needs to be added to the right- 
hand side of the equation to account for the external potential [Chandler 
(1987), Zwanzig (2001)]. Therefore, to describe the movement of a subdiffu- 
sive particle under the presence of an external potential U(x), we have the 
following model. 

The model for subdiffusion under a general potential U(x): 

m^p. = /* v { u )K H {t -u)du- U'(x(t)) + F H (t), 

Ctt J — oo 

x(t) = / v(s) ds, 





which is the companion of (2.7). 

The harmonic potential U(x) = ^rmfjx 2 , where m is the mass of the par- 
ticle and the constant ip reflects the potential's strength, is of particular 
importance because when the movement is confined to a short range, such 
as the movement within a protein as we shall see in Section 5, the underly- 
ing potential function can often be adequately approximated by a harmonic 
one. Under such a harmonic potential, the model is as follows. 

The model for subdiffusion under a harmonic potential: 

dvlt) /"* 
m—^- = -( / v(u)K H (t -u)du- mipx(t) + F H (t), 
dt J-oo 

(3-1) 

x(t) = / v (s) ds. 
Jo 

An equivalent expression in the more familiar probability notation is 
dx(t) = v(t)dt, 

(3.2) mdv(t) = -((J v{u)K H (t-u)du^jdt 



m 



i/)x(t) dt + v / 2(k B T dB H (t). 



Remark 4. Compared with the model (2.9) for the movement of a free 
particle, the above model (3.2) has two distinctive features. First, the pres- 
ence of the external potential changes the model from a single equation to 



MODELING SUBDIFFUSION WITHIN PROTEINS 



11 



a set of two coupled equations, and correspondingly, the solution to model 
(3.2) is a two-dimensional process. Second, for a free particle, since there 
is no external potential to bound its movement, the displacement process 
x(t) cannot be stationary [which is manifested in Theorem 2.2, where the 
variance Var(x(i)) oc t 2 ~ 2H — > oo, as t —> oo]. On the other hand, for model 
(3.2), since the harmonic potential always pulls the particle back to the ori- 
gin [because of the term — mipx(t)], the displacement process x(t) can be 
stationary, and thus, we can talk about the stationary mean and variance 
of x(t), as we shall see next. 

The next theorem, whose proof is deferred to the Appendix, gives the 
solution to model (3.2), which describes the subdiffusive motion under a 
harmonic potential. 

Theorem 3.1. Under the pathwise interpretation of dBjj(t) for 1/2 < 
H < 1, the solution to equation (3.2) is 



/oo 
p(t-u)dB H (t 
-oo 

/oo 
p'(t-u) dB H { 
-oo 



where the deterministic function p{t) is defined as 



OO 1 

e~ ituJ ^- dw . 

mip — muj 2 — iuJ^Kjj{uj) 



Furthermore, the solution (x(t),v(t)) is a stationary bivariate Gaussian pro- 
cess with mean function E[x(t)] = E[v(t)] = and covariance functions given 
by 

E[x(s)x{s + 1)} = E[x(0)x(t)] 



2-7T j-oo \mip — muj 2 — iu)C,KJj{u)\ 2 
E[x(s)v{s + t)} = E[x(0)v{t)} 

k ^ r e -»> , ^ 



27T j-oo \mtp — mu} 2 — iu}C,K^{ 



12 



E[v(s)x(s + <)] = E[v(Q)x(t)] 



k ^ r e -u» ^Mt ^ 



27r 7-oo \mtp — mw 2 — iu^Kjj{. 



12 



E[v(s)v(s + 1)] =E[v(0)v(t)] 
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where the expressions of Kh{u) and Kjj(uj) are given by (2.11) and (2.12), 
respectively. 

Remark 5. It is straightforward to verify that in the limiting case of 
tp — > 0, the result of E[v(0)v(t)] in the above theorem converges to the ex- 
pression of C v (t) in Theorem 2.1. This says that as the harmonic potential 
becomes weaker and weaker, the particle will behave more and more like 
a free particle, and in the limit the movement reduces to that of a free 
particle. 

One special case of model (3.1) is when the acceleration term mdv{t) / dt is 
negligible. This corresponds to the so-called overdamped condition in physics 
[Van Kampen (2001)], where the friction in the system is very large, causing 
the acceleration of the particle to be negligible. In the overdamped scenario, 
the acceleration term drops out, and the model changes to the following. 

The model for subdiffusion under a harmonic potential and the over- 
damped condition: 

(3.3) mipx(t) = -C f v{u)K H {t-u)du + F H (t), x{t) = [ v(s)ds. 

J-oo ' JO 

It can be rewritten in the more familiar probability notation as 
dx(t) = v(t)dt, 

(3-4) 

mipx(t) dt = -( (^J v{u)K H {t - u) du^j dt + ^2(k B T dB H (t). 

The next theorem (with proof given in the Appendix) solves equation 

(3.4) . We will use the solution to explain the experimental data in Section 5 
because in biological systems, such as within a protein, the friction is usually 
large and the overdamped condition usually holds. 

Theorem 3.2. Under the pathwise interpretation of the stochastic inte- 
gral for 1/2 < H < 1, the solution to equation (3.4) is a stationary Gaussian 
process with mean E[x(t)\ =0 and covariance function 

(3.5) a x (t) = E[x(0)x(t)] = ^I E ^ 2H (-(t/r) 2 - 2H ), 

imp 

where the constant 

1/(2-2//) 

(3.6) 



/ Cr(2ff + i) y 
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and E a (z) is the Mittag-Leffler function [see Erdelyi et al. (1953), Chapter 
18] defined by 

oo 

E a {z) = ^z k /T(ak + 1). 

k=0 

Remark 6. The Mittag-Leffler function generalizes the exponential func- 
tion in a natural way. When H — > 1/2, the Mittag-Leffler function in (3.5) 
reduces to the exponential function, and a x (t) = (kBT/mtp) exp(—(mtp/()t), 
recovering the classical Brownian diffusion result. 

4. Physical basis of the model. We shall apply the results in the previous 
two sections to explain the nanoscale subdiffusive motion observed within 
proteins. But before doing so, we will study in this section the physical 
foundation of the model, since a key requirement for biophysical models is 
that, in addition to satisfying fundamental physical laws, they must have a 
sound physical basis. 

4.1. The thermal dynamic requirement for a free particle. The law of 
thermal dynamics [Chandler (1987), Hill (1986), Reif (1965)] requires that, 
for a free particle, the (equilibrium) stationary variance of its velocity should 
be ksT/m, where m is the mass of the particle. The next theorem (whose 
proof is deferred to the Appendix) verifies that indeed our model (2.9) for 
the free particle satisfies this thermal dynamic requirement. 

Theorem 4.1. Under model (2.9), the stationary variance of the veloc- 
ity Var[f (0)] satisfies 

Var[«(0)] = — for all 1/2 <H<1. 

4.2. The thermal dynamic requirement for the movement of a particle 
under harmonic potential. For particles moving under a harmonic potential 
U(x) = ^mifjx 2 , the law of thermal dynamics asserts that the equilibrium 
(stationary) variance of the displacement should be The next theorem 
(whose proof is given in the Appendix) confirms that our model (3.2) for 
harmonic potential indeed satisfies the thermal dynamic requirement. 

Theorem 4.2. Under model (3.2), the stationary variance of the dis- 
placement Var[x(0)] satisfies the thermal dynamic requirement of 



Var[x(0)l = for all 1/2 < H < 1. 
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For model (3.4), which describes subdiffusion under harmonic potential 
and the overdamped condition, the results of Theorem 3.2 imply 

(4.1) ^(0) = Var[x(0)] = ^E 2 _ 2H (0) = for all 1/2 < H < 1. 

mip mip 

We thus have the following: 

Theorem 4.3. Under model (3.4), the stationary variance of the dis- 
placement Var[x(0)] satisfies the thermal dynamic requirement of 

Var[z(0)l = -2-r for all 1/2 < H < 1. 

mip 

4.3. Deriving the model from a system of interacting particles. In this 
subsection we will demonstrate that the model in Section 3 can be derived 
from the physical microscopic interaction between the particle under study 
and its surrounding media; in particular, the model can be derived from a 
Hamiltonian system consisting of the particle and its surroundings. For more 
general discussion about the Hamiltonian and GLE, see Zwanzig (2001). 

We start the derivation from the Hamiltonian [Corben and Stehle (1995)] 
of the particle 

p 2 1 

H s = - h -mipx 2 , 

2m 2 

where p = mv is the momentum, and x is the displacement. Here the Hamil- 
tonian, which is essentially the total energy of the particle, consists of the 
kinetic energy p 2 /(2m) = mv 2 /2 and the potential energy, which is mipx 2 /2 
under the harmonic case. The surrounding media, consisting of N small 
molecules, has its own Hamiltonian (total energy) 

N ( v 2 1 / v x 2 



where N, on the order of 10 23 , is the total number of molecules in the media, 
?7ifo is the (common) mass of an individual molecule in the media, pj and 
qj are respectively the momentum and location of the jth molecule, Uj is 
the oscillation frequency of the jth molecule (as each individual molecule 
oscillates in the media), and the term (qj —^fjx/ui 2 ) captures the interac- 
tion between the particle of interest and the individual molecules in the 
media, where jj is the interacting strength between the particle and the jth 
molecule. 

The total Hamiltonian (energy) of the entire system (i.e., the particle plus 
the media) is thus H s + Hb- Once the total Hamiltonian is given, the classical 
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theory of mechanics [Corben and Stehle (1995)] states that the motion of the 
particle, as well as that of the individual molecules, is given by 



(4.3) 
(4.4) 



dx = d(H s + H B ) dp = d{H s + H B ) 
dt dp dt dx 

dq 3 _ d{H s + H B ) dp, _ d(H s + H B ) 



dt 



d Pj 



dt 



dqj 



which is a set of coupled differential equations. The exact expressions of H s 
and Hb reduce (4.3) and (4.4) to 



(4.5) 
(4.6) 



dx 


- IL 


dp 


~dt 


m' 


It 


dqj_ 


= — 


dpj_ 


dt 


m 6 ' 


dt 



N 



-mtpx + ^2jjm b (qj 



3=1 



-rribWjqj + m b -fjx, 



2 X I ' 



from which we can first express qj and pj in terms of x(t) and their initial 
values: 



qj(t) = qj(0) cos(uJjt) + t ' J — sm(u>jt) H — - / x(s) s'm(ujj(t — s)) ds. 

rribUJj ujj Jo 



Applying an integration by parts on it gives 



qj(t) - ^x(i) 



^■(O)-^x(O) 



/ \ Pj'(O) • / N 

cos[Ujt) H — sin(w,-i) 

TTlifLOj 



— -^L cos(ujj(t — s)) dx(s). 
u>]Jo 

Taking this expression into (4.5), we obtain 

d?x f t 
m—7r = —mipx(t)— / Jit — s) dxis) + G(t), 
dt z Jo 



where 



N 7 2 

J(t) = nib -^2 cos(ujjt)i 

3=1 ^3 

N 



G(t) = m blj ( 9j(0) - ^Jx(O) ) cos(wjt) + J] 77Pi(0) sin^i). 



^7 



N 



1 W 7 
J=l J 



Upon making the association 

Cocm b , (K(t) = J(t), 
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we reach the GLE model 

~dt 



m—p- = -mipx(t) - C / Kit- s)v(s)ds + G(t), x(t) = [ v(s)ds. 

Jo Jo 



Remark 7. If we let the harmonic potential rmpx 2 /2 become weaker 
and weaker, then in the limit of ip — > 0, the model (2.5) of a free particle is 
recovered. 

The randomness of G(t) comes from the fact that the initial values of q = 
(<Zi(0)i 92(0), . . . , gjv(0)) and p = (pi(0), . . . ,pjv(0)) have a thermal dynamic 
distribution 

/(p,q) ocexp^-^0, 

which leads to 

E[p 3 (0)] = E[ qj (0) - 7;z(0)/^ 2 ] = for all j, 



^[Pj (0)] = "i 6 A;bT, S[(«(0) - 7i»(0)/ 



7 V. w /J '•"0" J rs-^ 1 WHiy^ ) 13 \ / 1 7 / J 2' 

^bj - 7^(°) M 2 )] = for a11 * and i) 

implying 

N 7 2 

£[G(t)G(«)] = k B Tm b y2^ cos(ujj(t - s)) = k B TCK(t - s), 

0=1 J 

which exactly recovers the fluctuation-dissipation relationship (2.6). 

The underpinning of a fractional Gaussian memory kernel. Because N, 
the total number of molecules in the media, is large, the memory kernel 
K(t) oc ■^J2f=i7j cos ( LJ ji)/ u ''j i s essentially given by 

1 roc E[-i 2 \uj] 



K{t) oc E 



y 

— =■ cos(ujt) 



w 



■ cos(ut)g(u;) duj, 



where g(-) is the probability density function of the molecules' oscillation 
frequencies. If 7 and g(u>) are such that 

(4.7) 3Z£W„' 



,1~2H 

then 

|2fl-2 



K(t)oc\t\^ H ~ z otK H (t), 



MODELING SUBDIFFUSION WITHIN PROTEINS 



17 



giving rise to the memory kernel (2.10) of the fractional Gaussian noise. 
Many scenarios can lead to equation (4.7). For example, if the interacting 
strength 7 is a (deterministic) power function of w. 7 2 ~ oj 3 ~ 2H , and the 
distribution g is roughly uniform over the spectrum, then (4.7) would hold 
approximately. Another example is 7 and u being independent and g having 
an (approximate) power tail g(u>) ~ u> 3 ~ 2H exp(— au>) with a very small a; 
then (4.7) would also hold approximately. 

5. From theory to experiments. A recent single-molecule experiment 
[Yang et al. (2003)] studied a protein-enzyme compound Fre, which is in- 
volved in the DNA synthesis of E. Coli. As shown in Figure 1, Fre contains 
two subunits: FAD and Tyr. Because the 3D conformation of Fre sponta- 
neously fluctuates over time, the (edge-to-edge) distance between FAD and 
Tyr varies. This distance is one dimensional; its fluctuation provides infor- 
mation about the conformational dynamics of Fre. To experimentally probe 
this one-dimensional distance fluctuation, Fre is placed under a laser beam. 
The laser excites FAD to be fluorescent. By recording the fluorescence life- 
time of FAD, one can trace the distance between FAD and Tyr, because at 
any time t the fluorescence lifetime X(t) of FAD is a function of the one- 
dimensional distance [see Gray and Winkler (1996), Moser et al. (1992)] 

(5.1) \(t) = k e? ( - x <* +x ®\ 

where k$ and (3 are known constants [Moser et al. (1992)], x eq is the mean 
distance, and x(t) with mean is the distance fluctuation at time t. 

To model x(t), we first note that the external potential experienced by 
the fluctuating subunits is well approximated by a harmonic one, U(x) = 
mipx 2 /2, because the movement is confined within the short range of Fre. 
We shall see in Section 5.3 that this close approximation is well tested in 
the experiment. 

5.1. Testing the autocorrelation structure of the model. With the har- 
monic potential, people used to model x(t) as a Brownian diffusion process 
m4fv(t) = —(v(t) — mipx{t) + F(t),x(t) = Jqv(s)cIs, or by its overdamped 
version mipx(t) = —£v(t) + F(t),x(t) = J^v(s)ds, where F(t) is the white 
noise. 

The nanoscale single- molecule experimental data of A(i), unlike the tra- 
ditional population experiments, provides the means to test the model. One 
can calculate the empirical autocorrelation function of X(t) from the exper- 
imental data and compare it with the theoretical autocorrelation function 
from the model. The autocorrelation function is used as the test statistic be- 
cause the experimentally recorded fluorescence lifetime is actually the true 
X(t) plus background and equipment noise. Doing an autocorrelation effec- 
tively removes the noise (since the noise is uncor related). For a stationary 
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Fig. 2. Autocorrelation function of the fluorescence lifetime X(t). The open circles repre- 
sent the empirical autocorrelation calculated from the experimental data. The dashed line is 
the best fit from the classical Brownian diffusion model. The solid line is the fit [H — 0.74, 
£/(mip) =0.40, (3 2 kBT/(mip) = 0.81/ from our model (3.4), agreeing well with the data. 

Gaussian process x(t), it is straightforward to calculate the autocorrelation 
function of A(t) from equation (5.1), 

(5.2) cov(A(0),A(i)) = ^e^+^^toif/ftW _ i) ; 

where C x (t) = cov(x(0), x(t)). Figure 2 shows the empirical autocorrelation 
function (the open circles) compared with the best (least-square) fitting 
from the Brownian diffusion model (the dashed curve). A clear discrepancy 
is seen. 

The solid line in Figure 2 is the result from modeling x(t) by the subd- 
iffusive process (3.4) under the harmonic potential. The curve is fitted by 
using the Hurst parameter H = 0.74, expression (5.2) and the result of Theo- 
rem 3.2 [with o~ x {t) replacing C x {t) in (5.2)]. A very close agreement with the 
experimental autocorrelation function is seen. Here the overdamped model 
(3.4) is applied to explain the data, since the movement within a protein is 
subject to the overdamped regime. 

5.2. Testing higher-order correlation functions. To check our model, we 
make predictions about the distance fluctuation and test whether these pre- 
dictions can be confirmed by the experiments. The first set of predictions in- 
volves higher-order autocorrelation functions because they are very sensitive 
to distinguishing models [Mukamel (1995)]. With the values of the fitting pa- 
rameters fixed to those in Figure 2, we compute from the model the predicted 
three-step and four-step autocorrelation functions i?[AA(0)AA(ii)AA(ii + 
t 2 )} and £[AA(0)AA(ii)AA(ii + t 2 )AA(fi + t 2 + t 3 )), where AA(t) = X(t) - 




Fig. 3. Higher-order autocorrelation functions ofX(t). (a), (b) and (c): The experimen- 
tally obtained autocorrelation functions E[AX(0)AX(t)AX(2t)], £7[AA(0)AA(2t)AA(3i)] 
and E[AX(0)AX(t)AX(2t)AX(3t)] overlaid with the model predictions for various t. The 
theoretical curves from the model (3.4) are calculated using the same parameter values as 
in Figure 2 [H = 0.74, C/(mi>) = 0.40, and /3 2 k B T / (rmp) = 0.81] . 



E[X(t)], and compare them with their experimental counterparts. The exact 
expressions for the three-step and four-step autocorrelation functions are 
given in the Appendix. 

Figure 3(a) shows the evenly spaced three-step autocorrelation function 
E[AX(0)AX(i)AX(2t)] as a function of time t; Figure 3(b) shows the un- 
evenly spaced three-step autocorrelation function i?[AA(0)AA(2i)AA(3i)] 
as a function of time t; Figure 3(c) shows the evenly spaced four-step au- 
tocorrelation function £'[AA(0)AA(t)AA(2t)AA(3t)] as a function of t. The 
theoretical curves (the solid lines) in Figure 3 are calculated from model 
(3.4) using the parameter values obtained from the fitting in Figure 2. In all 
cases, the close agreement between the theoretical curves (the solid lines) 
and the experimental values (the open circles) is seen. 

The second prediction from the model is time-symmetry. For any t\ and 
t 2 , the model predicts E[A\(0)AX(t 1 )AX(t 1 + t 2 )] = E[AX{0)AX(t 2 )A\(ti + 
£2)]; which can be straightforwardly seen from the formulas in the Appendix. 
It says that if our model is true, then one can swap the order of the time lags 
without changing the correlation value. This can be tested by taking ti = 
t, t 2 = 2t and plotting the experimentally obtained three-step correlation 
E[AX(0)AX(t)AX(3t)\ against £[AA(0)AA(2f)AA(3t)] for various t. A 45° 
line is predicted by the model. The experimental plot in Figure 4 indeed 
confirms the prediction. 

5.3. Testing the harmonic potential. As a final check of our model, we ask 
if the two important model assumptions of harmonic potential and fractional 
Gaussian memory kernel (2.10) can be directly verified from the experiment. 
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0-00 0.01 0.02 0.03 

£'|AX(0) AX(2o AX(3f)] 

Fig. 4. j4 test for time-symmetry. The experimental three-step correlations 
E [AA (0) A A(t) AA (3t)] and £[AA(0)AA(2t)AA(3f)] plotted again each other for 
various t. A 45° line is predicted by our model. 



Another recent single-molecule experiment [Min et al. (2005)] indeed con- 
firmed the assumptions. This experiment studied a protein complex formed 
by fluorescein (FL) and monoclonal antifluorescein (anti-FL). See Figure 5. 
Similar to Fre, this complex contains two substructures Tyr and FL, between 
which the distance fluctuates over time. Using exactly the same experimen- 
tal technique as in the previous Fre experiment, the distance fluctuation can 
be probed from the fluorescence lifetime of FL (upon placing the complex 
under a laser beam). 

This latter experiment has identical settings as the previous one, ex- 
cept that it has much higher signal-to-noise ratio, and thus provides higher 
resolution data that allows the model assumptions to be further tested. 
The high resolution data on X(t) first enables x{t) to be reconstructed 
from (5.1) through a local binning (kernel) average. See Figure 6(a). The 




Fig. 5. The crystal structure of the FL and anti-FL complex. The two substructures Tyr 
and FL are highlighted. This new experiment has idential settings as the previous Fre ex- 
periment (Figure 1), but it has much higher signal-to-noise ratio that enables experimental 
testing of the key assumptions of our model. 
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Fig. 6. (a) The distance fluctuation x(t) reconstructed from (5.1), and the corre- 
sponding empirical distribution P(x). (b) The estimated empirical potential function 
U(x) — —kBT\og(P(x)), the open circles, compared with the harmonic potential, the 
dashed line. 



empirical equilibrium (stationary) distribution P{x) of x(t) is then ob- 
tained from the histogram of all the x{t). According to thermal dynam- 
ics, the equilibrium distribution P{x) and the potential function U(x) is 
linked by P(x) oc exp(— [/(^/(feeT)). A natural estimate for U(x) is thus 
U{x) = —kBTlog(P(x)), which is shown in Figure 6(b). A harmonic po- 
tential, the dashed line in Figure 6(b), is seen to fit U(x) very well, hence 
confirming the validity of our earlier assumption. 

5.4. Testing the fractional Gaussian memory kernel. To test the frac- 
tional Gaussian memory kernel (2.10), first we calculate from the recon- 
structed x(t) the experimental autocorrelation function C x (t) = cov(x(0),x(t)) 
and compare it with the theoretical expression (3.5) from Theorem 3.2. Fig- 
ure 7(a) shows the comparison, where the result from our model (the solid 
line) agrees well with the experimental values (the open circles). 

Furthermore, from the overdamped GLE mipx(t) = v(u)K(t — 

u) du + G(t),x(t) = Jqv(s)cIs with an arbitrary memory kernel K, we can 
establish a one-to-one correspondence between the Laplace transform of 
C x (t) — cov(x(0), x(t)) and the Laplace transform of the memory kernel 
K(t), as shown in the following theorem, which then allows us to recover 
K{t) from the experimental C x {t) = cov(x(0), x{t)) and compare it with the 
assumed fractional Gaussian noise memory kernel. 

Theorem 5.1. For the overdamped GLE mi/jx(t) = —( f^v^Kft — 
u)du + G(t),x(t) = jQv(s)ds with a general memory kernel K(t), there is a 
one-to-one correspondence between the Laplace transform C x (s) = 
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Fig. 7. (a) The autocorrelation functions of x(t) calculated from the experimental data, 
the open circles, compared with the Mittag-Leffler expression (3.5) from our model, the 
solid line, (b) The experimentally determined memory kernel compared with the memory 
kernel of the fractional Gaussian noise. 



J Q e C x (t)dt of C x (t) = cov(x(0), x(t)) and the Laplace transform &(s) 





Jq°° e st K(t) dt of the memory kernel K(t): 



(5.3) <t x (s) = kBTC ,* ( f> K^-^ m ^ x{s) 



mip mip + Cs&(s)' C IzbT — rmps<£ x (s) 

Therefore, knowing C x (t) allows the recovery of K(t) in the Laplace space. 



The proof of the theorem is given in the Appendix. With this theorem, 
using the empirical autocorrelation function C x (t), we can determine (in 
the Laplace space) the memory kernel experienced by the protein in the 
experiment. Figure 7(b) shows, in the Laplace space, the experimentally 
obtained memory kernel from (5.3) compared with the Laplace transform of 
the fractional Gaussian memory kernel &h(s) = Jo° e~ st KH(t) dt = T(2H + 
l)^ 1 " 2 ^, which is a power law. Close agreement is seen, which verifies the 
second key assumption of our model. 



6. Discussion. To explain the experimentally observed subdiffusion phe- 
nomenon, we formulate in this article a stochastic model by incorporating 
fractional Gaussian noise into the generalized Langevin equation framework. 
The resulting stochastic integro-differential equations driven by fractional 
Brownian motion are nonstandard. We study in detail these model equa- 
tions. Using the analytical results, we show that the model leads to a satis- 
factory account for subdiffusion. 
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The model, in addition, has three attractive features: (1) The under- 
lying theory is simple. First, compared with the classical Brownian diffu- 
sion theory, the model has only one more parameter: the Hurst parameter 
H . Second, the model offers analytical tractability. For instance, under the 
harmonic potential, closed form expressions of the displacement covariance 
function are obtained in Theorems 3.1 and 3.2. (2) The model, derivable 
from a Hamiltonian consideration, has a sound physical basis, which is an 
important requirement for biophysical models. (3) The theoretical results 
from the model agree well with the experimental data. Not only are the 
model predictions confirmed by the experiments, but also each key model 
assumption is directly verified in the experiments. 

The successful application of our model to explain subdiffusion only ex- 
emplifies one instance of the numerous and growing stochastic modeling 
opportunities in nanoscale biophysics. Many interesting problems remain to 
be explored. 

1. Existence, uniqueness and solution under a general potential. In this 
paper we solved the model equations (GLE with fGn) for the harmonic po- 
tential U(x) = rmpx 2 /2. For a general potential U(x), the equation becomes 



A natural follow-up question is as follows: under what conditions does there 
exist such a bivariate process (x(t),v(t)), and when is it unique, for example, 
in the weak sense? Furthermore, if such a bivariate process exists and is 
unique, then how might one solve it, at least numerically? The answers to 
these questions are directly related to many biological and chemical systems, 
since many such biophysical and biochemical processes are subject to general 
potentials. 

2. First passage time and rare events calculation. Many biological events 
are associated with the first passage time problem; for example, the comple- 
tion of an enzymatic reaction corresponds to the first time that the enzy- 
matic system escapes an energy barrier [Risken (1989), Van Kampen (2001)]. 
Consequently, calculating the distribution of the first passage time is of im- 
mediate biological applicability. For instance, the distribution of the first 
time that x(t) reaches a given level from equation (6.1) directly relates to 
the understanding of enzymatic reactions in a sluggish protein environment 
[Min et al. (2005)]. Furthermore, these first passage time problems in biology 
usually correspond to the crossing of a very high barrier. An interesting open 
problem is thus to investigate how simulation, such as importance sampling, 



dx(t) 



v(t)dt, 





U'(x(t)) dt + ^Mk^TdB H {t). 
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or large deviation techniques can be applied here to assess or approximate 
the probability of high barrier crossing. 

3. Interacting particle systems. It is seen in Section 4 that our model can 
be derived from a system of interacting particles. This type of microscopic 
picture is not unusual for biological systems, as biological events, such as 
gene expression, tend to resolve from the interaction between many small 
units. It is thus interesting to see how the modern understanding of interact- 
ing particle systems can be extended to biological systems. For our model, 
in particular, a linear coupling between qj and x in the Hamiltonian (4.2) is 
assumed in the derivation. It is desirable to relax this assumption to extend 
the microscopic derivation to more general interaction terms. 

The booming field of nanoscale (single-molecule) biophysics has attracted 
much attention from biologists, chemists and physicists, as it projects a 
bright picture for new scientific discoveries. It also presents many interest- 
ing and challenging problems for stochastic modelers due to the stochastic 
nature of the nanometer world. It is our hope that this article will generate 
further interest in applying modern probabilistic and statistical methodol- 
ogy to interesting biophysical and scientific problems. 

APPENDIX: PROOFS AND DERIVATIONS 



Proof of Theorem 2.1. The pathwise interpretation oidB H {t) allows 
us to apply the technique of analyzing classical integro-differential equations 
to solve (2.9). Suppose v(t) is the solution. Then applying a Fourier trans- 
form to both sides of equation (2.9), we know that v(u>) = f^°^e %tu v(t)dt 
must satisfy 

/OO 
e itu} dB H (t). 
-OO 

The unique solution of the above equation is 

/OO 
e ltuJ dB H (t)/[(K+(u;) - imuj}. 
-OO 

An inverse Fourier transform gives 

1 /*oo /*oo 

v(t) = — e~ ituJ v(uj) du = ^2(k B T / r(t - u) dB H (u). 

^7T J — oo J -co 

Since r(t) is deterministic, it is straightforward to verify that v(t) given 
above is indeed a finite Gaussian process with zero mean. To calculate the co- 
variance function, we use the fact [cf. Duncan, Hu and Pasik-Duncan (2000)] 
that, for H > 1/2 and deterministic functions / and g, 

E 



f{u)dB H {u)- J g{u)dB H {u) 

H{2H - l)\u - v\ 2H - 2 f{u)g{v) dudv, 
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which provides 

E[v(t)v(s)} 

/oo poo 
/ H(2H - l)\u - v\ 2H ~ 2 r(t - u)r(s - v) dudv 
-oo J —oo 

/•oo roo 

= k B TC, I I r(t — u)r(s — v)Kh(u — v) dudv. 



oo J — oo 



A change of variable of y = t — u, z = s — v gives 

r-OO POO 



E[v(t)v(s)] = k B T( 



r{y)r{z)K H {t - s -y + z)dydz 



oo J —oo 



= E[v(0)v(t - s)]. 

Therefore, v(t) is a stationary Gaussian process. To obtain the stationary 
covariance C v (t) = E[v(0)v(t)], we compute C v (uj) = e ttu; 'C v (t) dt. 

From (A.l), it follows by a change of variable of y = t — u, z = u — v that 



e ltw C v (t)dt 



/OO POO POO 

dte itu) / r{t-u)r{-v)K H {u-v)dudv 
-oo J~ oo J — oo 



oo 
oo 



k B TQ[ / r(y)e tyw dy ) [ / r(-v)e tvu dv 



K H (z)e izuJ dz 



By the definition of r(t), the above expression is simplified to 

C v (cj) = k B T(KH(uj)/[((K+(Lj) - ium)((K+(-uj) + ium)], 
which for u G R can be further simplified to 
(A.2) C v (lo) = k B TtK H {uj)/\tK+{uj) - iujm\ 2 

because K^(—u>) is the complex conjugate of K^(uj) as the kernel function 
K}j{t) is real. □ 

Proof of Theorem 2.2. Since E[x(t)} = Jq E[v(s)]ds = 0, it follows 
that, for t > 0, 



Var[x(t)] =E[x 2 {t)} 
The right-hand side is 
2 



t rt 



o J o 



E[v (s)v (u)) duds = 2 / / C v (u) duds. 
Jo Jo 



t 

o Jo 



C v {u) duds = 2 
= 2 



t rs 



1 

o [2tt 

t MM 



o Jo 



vr Jo 



e- tuu C v (Lo)duj 



cos(ulu)C v (uj) duj 



duds 
duds. 
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Applying Fubini's theorem twice simplifies it to 

pt rs 2 r°° 1 

2 / / C v (u) duds = — / — zC v (u))(l — costuj) dw. 
Jo Jo n Jo ^ 

Plugging in the result of Theorem 2.1, we obtain 
2 mi _ 2 f° 1 k B TCK H {u) 



7T 



/o W I QKjjiuj) — imuj\ z 



2 /*°° 

= - / (2fc B rcr(2F + l)sin(i?7r)a;~ 1 ~ 2// (l-costa;)) 

x (m 2 cj 2 + 2r(2if + 1) cos(HTr)m(uJ 2 ~ 2H 

+ C 2 T 2 {2H+l)uj 2 - m y l du. 

A change of variable i] = tu> gives 
E[x 2 (t)] 

= - t 2 ~ 2H / (2k B TQT{2H + 1) sin(i2'7r)7/- 1 - 2ff (1 - cost/)) 
vr Jo 

x {ra 2 T] 2 t~ m + 2r(2# + 1) cos(H-K)mCr] 2 ~ 2H t~ 2H 

+ ( 2 T 2 (2H+l)r 1 2 - 4H r 1 dri. 

As i — ► oo, by the dominated convergence theorem, the above integral con- 
verges to 

2k B Tsm(H-K)r]- l ~ 2H (l - cost/) /c b T sin(2iT7r) 



/ 

Jo 



(F(2H + l)r] 2 - 4H 1 ( 2H(l-2H)(2-2H)' 

Therefore, as t — ► oo, the mean-squared displacement 

2-2H k B T sm(2Hir) 



E[x 2 {t)}/e 



C ttH{\-2H){2-2H)' □ 



Proof of Theorem 3.1. Applying a Fourier transform to equation 
(3.2), we know that x{u) = Jf^e'^ 'x(t) dt and v{u) = J^^vtf) dt must 
satisfy 

v(uj) = —iujx(Lo), 



-imujv{uj) = -(v(lj)K+(uj) - rmpx(u) + y / 2(k B T / e ituJ dB H (t), 

J — oo 

which has the unique solution 

x(uj) = — , - + / e dB H (t), 

v{u) = — ~ 2 — . ^ +( : / e dB H (t). 
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An inverse Fourier transform gives 

/oo 
p(t-u)dB H (u), 
-oo 



1 — oo 
roo 



/oo 
p'{t-u) dB H {u). 
-oo 



With p(t) being deterministic, it is straightforward to verify that (x(t),v(t)) 
given above is a finite stationary bivariate Gaussian process with zero mean. 
For the covariance function E[x(Q)x(t)], we have 

/oo roo 
/ H(2H - l)\u - v\ 2H ~ 2 p(t - u)p{-v) dudv 
-oo J — oo 

/OO /"OO 
/ p(t — u)p{—v)Kn(u — v) dudv. 
-oo J — oo 

The Fourier transform of the above equation is 

COO 

e itu} E[x{0)x(t)]dt 



■oo 



k B T( I dte ltuJ / p(t-u)p(-v)K H (u-v)dudv 

-oo J —oo J —oo 



k B T( 



oo 



y = t — u, z = u — v. 
By the definition of the above expression is simplified to 

e ituj E[x(0)x(t)} dt = k B T(K H (uj)/\rmp - mto 2 - iu<i^(w)| 2 . 

oo 



oo 



Thus, 



E[x(0)x(t)] = r e- itw K H {u)l\rml) - muj 2 - iu(K+{u)\ 2 du. 

27T J —oo 

The expressions of E[x(0)v(t)], E[v(0)x(t)] and E[v(0)v(t)] can be obtained 
similarly. □ 

Proof of Theorem 3.2. Following the proofs of Theorems 2.1 and 3.1, 
applying the Fourier method on equation (3.4) and some detailed calcula- 
tions afterward yield 



/oo 
p{t-u) dB H {u 
-oo 

/oo 
p!{t-u) dB H {i 
-oo 
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where 



1 



oo 



H(t) = — I e~ it0J dw, 

2ir J-oo mtp — iu(,Kfj{uj) 

from which we know that E[x(t)] = 0, and that the Fourier transform a x {uo) 
IZo eitU(7 x{t) dt = X!^ e itui E[x(<S)x{t)} dt satisfies 

a x (u) = k B T(K H (oj)/\mij; - iuj(K^(uj)\ 2 . 
Using the expressions (2.11) and (2.12), we have, for to > 0, 
a x (u) = (2k B T(T(2H + 1) sin(i7vr)tj 1 " 2H ) 

x (m 2 tp 2 - 2mip(T(2H + 1) cos(Htt)uj 2 ~ 2H 

(A.3) 

+ C 2 r 2 (2if + i) w 4 " 4 ^)- 1 

_ k B T 2sm(Hir)(TLv) 2 - 2H /u; 



rrnp 1-2cos{Htt)(tuj) 2 - 2H + (™) 4 - 4ff ' 

where the last equality uses the definition of r in (3.6). 

To establish (3.5), we only need to show that the Fourier transform 
of ^-E 2 -2H(-{t/T) 2 ~ 2H ) is exactly equal to a x (u), that is, 
^ 00 e- itw {k B T/m^)E 2 -2H{-{t/T) 2 ~ 2H ) dt = a x (u), which by (A.3) reduces 
to show that, for to > 0, 

roc 

2 / COs(tL0)E2-2H(-(t/T) 2 - 2H )dt 

(A.4) J ° 

2sm(Hir)(TLj) 2 ~ 2H fu 



1 - 2cos(ii^)(™) 2 - 2 " + (rw) 4 - 4 ^ ' 

The Laplace transform of the Mittag-Leffier function has been given in 
Erdelyi et al. (1953), Chapter 18 as 



OO 1 

e pt E a {-{t/r) a )dt-- 



P 1 + (rp)' 
Taking p = ito in the above formula gives 

2 cos(tuj)E 2 _ 2H (-(t/T) 2 - 2H )dt 
Jo 



2 Re 



1 1 



i0J 1 + {iTu)~£- 2H ) 



2 /l 1 
-Re 



LO 



i 1 + ( rw )-(2-2ff) e -i(l-tf)* 



2 ( 1 
-Re 



LO 



i + (tlu)-( 2 - 2H ) (-t cos(fl7r) + sin(tf vr)) 
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_ 2 sm(HiT)(TUj) 2 - 2H 
"wl- 2cos(^7r)(Tw) 2 - 2 ^ + (rw) 4 " 45 ' 

which is exactly (A. 4). The proof is thus completed. □ 

PROOF of Theorem 4.1. From the result of Theorem 2.1, we have 
Var{v{0)]=E[v 2 (0)] 

1 poo 

= — / k B T(K H (u)/\(K+(uj) - imuj\ 2 doj 

2TT J-oo 



oc 

oo 



= - / {2k B T(T(2H + l)sm(HTr)uj 1 ~ 2H ) 
vr Jo 

x(( 2 r 2 (2ff+l) W 2 - 4F + my 

+ 2m(T(2H + l)u; 2 ~ 2H cos(Ftt))" 1 dw. 
A change of variable r] = lo 2H gives 

VarW0n = ^ r Cr(2g+l)sin(^) dv 

[ { irH Jo mV + 2mCr(2iJ + l)cos(iJ7r)77 + C 2 r 2 (2F + l) 

Using the general formula f n °° ■> , » rfx , , 3 = — t— r , the above expression 

00 JO x z -\-2xy cos q>+y z ysin</>' ^ 

is simplified to 

Var[u(0)] = — for all 1/2<#<1, 
m 

which agrees with the thermal dynamic requirement. □ 

PROOF of Theorem 4.2. Theorem 3.1 implies 
Var[x(0)] 

= E[x 2 (0)} 
k B TC, f°° 
2ir J-oo 
k B T roc 



vr Jo 



poo 

/ (2Cr(2F + l)sin(F7r)w 1 - 2 ^) 
Jo 

x ([mip - muj 2 - (T(2H + l)u 2 ~ 2H cos{Htt)] 2 

+ [(T{2H + 1)oj 2 ~ 2H sin(i77r)] 2 )~ 1 du>. 



A change of variable 77 = rw, where r = [^r(2i7 + l)] 1 ' ' 2 2 ^ , gives 

Var[x(0)] = — r \2sm(HTr)r ] 1 - 2H ) 
vr Jo 
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(A.5) x ([mif) - mrf/r 2 - rf- 2H cos(Hir)] 2 

+ { V 2 ~ 2H sm(HTr)} 2 y 1 dr ] . 

Consider the complex valued function 

f(z) = [z(mi/> - mz 2 /r 2 - z 2 - 2H e iHn )}-\ 

It is straightforward to verify that it is analytic on the region defined by the 
boundary curve 

C = [1/R, R] U {Re 10 : < 6 < vr} U [-R, -1/R] U {e ie /R :0<6<tt}, 

where the real number R > 1, and [1/R, R] is the real interval between 1/R 
and R. 

It follows that J c f(z) dz = 0. But 



where 



[ f(z)dz = I + II + III + IV, 
Jc 



1= / [x{rmp-mx 2 /T 2 -x 2 - 2tl e ltl7I )}- l dx, 
Jl/R 

II = r[Re w {m^ - mR 2 e 2W /r 2 - #2-2/^(2-2//) ^tf 
Jo 



10 
r-l/R 

111= / [x(mip-mx 2 /T 2 -x 2 ~ 2H e tH7t )]~ l dx, 

J-R 



— mw — m — — - rr e 

R V R 2 t 2 R 2 ~ 2H 



-1 je 



IV-- 

We thus have I + III = - (II + 77). We can simplify I + III, II and 7Y 



as 



I + IH= j R \x(mi\> - mx 2 /t 2 - x 2 ~ 2H e lU *)Y x dx 
Jl/R 

rR 

„2 1 2 2-2H -iH-K\i-l 



+ / \x(m^-mx l /T l -x l - ltl e- lti *)Y v dx 
Jl/R 

R 1 2ism(HTr)x 2 ~ 2H 

dx 

1/r x [mip — mx 2 /r 2 — x 2 ~ 2H cos(i?7r)] 2 + [x 2 ~ 2H s'm(HiT)] 2 

II = i \\mi> - mR 2 e 2W /r 2 - ^-2/^(2-2//) jHvyi ^ 
Jo 

J2i0 A6(2-2H) \ -1 
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It therefore follows that 

rR ! 



2 s\-n{Hit)x 



2-2H 



1/r x [mip — mx 2 /r 2 — x 2 ~ 2H cos{Hn)] 2 + [x 2 ~ 2H sm(Hir)] 2 



dx 



mip — m- 



D 2id i9(2-2H) \ -1 



dO 



' R 2 T 2 R 2-2H 

- r{rn^-mR 2 e 2ie lT 2 -R 2 - 2H ^ 2 - 2H )e iH ^)- 1 de. 
Jo 

Letting R — > +oo provides (by the dominated convergence theorem) 

c 2iQ A9(2-2H) \ -1 

e e JHtx \ 



lim 

R— >+oo Jo 



mip — m 



R 2 T 2 R 2-2H 



de 



mip ' 



lim 

R— >+oo Jo 



(mip - mR 2 e 2ie /r 2 - #2-2*^(2-2*0^-1 dQ = ^ 



yielding 



1 



2sin(f/7r)x 



2-2_ff 



dx 



7T 



Iq x [mip — mx 2 /t 2 — x 2 2H cos(Hir)] 2 + [x 2 2H sin(ii~7r)] 2 mip 
Plugging this expression into (A.5), we finally obtain 



E[x 2 (0)] 



k B T 

mip 



for all l/2<iT<l. 



□ 



Proof of Theorem 5.1. Consider the function 

k B T( 



h(s) 



mip mip + (sK(s) 



where &(s) is the Laplace transform of the memory kernel K{t). The inverse 
Laplace transform h(t) of h(s) is given by Doetsch (1974), pages 4 and 148, 

1 

2ni 
1 

2^ 
2 

7T Jo 



h(t) = ^- I e st h{s)ds 

e lult h(iu>) du), oj = is 

00 

cos(tji) Re[h(iuj)] doj. 



Since 



Re[MW)] = ^^Re 



£(iu;) 



+ i(ujfi.(iuj) 



k B T(Re[8.{im)] 
\mip — i(u>&(— iu>)\ 2 
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and 

POO 

K(lo)= / e ituJ K{t)dt = 2Re[£(iw)], 



roo 

K+(u) = / e itu K(t)dt = A(-iu), 
Jo 



we obtain 



1 f°° 

h(t) = - / cos(u;t)A; B rCK(u;)/|m^ -<^ + (^)| 2 da;. 
vr Jo 

On the other hand, from the proof of Theorem 3.2, we know that for a 
general memory kernel K(t), the covariance C x (t) = cov(x(0), x(t)) is given 
by 

1 f°° 

C x (t) = — / e- ltu) k B T(K(uj)/\m4)-iuj(K + (uj)\ 2 duj 



1 

7T JO 



oo 

oo 



cos(o;t)A;BTCK(u;)/|m'i/' — iQuoK (o;)| cL>, 



which is identical to /i(t). Therefore, the Laplace transform <£ x (s) of C x {t) 
is /j(s), namely, 



Exact expressions of the higher-order autocorrelation functions of X(t). 
Since the fluorescence lifetime A(t) and the distance fluctuation x(t) are 
linked by 

A(t) = A;oe^ +x W), 

to calculate the higher-order autocorrelations of A(i), the following expres- 
sion is very useful. 

A useful expression. Suppose x(t) is a stationary Gaussian process with 
mean 0, and covariance function C x (t) = Cov(x(0),x(t)). Then for any t\ < 
h < ••• < t n , the expectation of E{e Ax ^ e Ax ^ ■ ■ ■ e Ax ^} , where A is a 
constant, is given by 

(A.6) E{e Ax ^e Ax ^ ■ ■ ■ e Ax ^} = expl -A 2 C x (0) + A 2 J2 C x {tj - U)\ ■ 

{ i<j J 

Proof. Since x{t\) + • — \- x(t n ) is Gaussian with mean and variance 
nC x (0) + 2J2 i<j C x (t j -U), it follows that e Ax ^ e Ax ^ ■ ■ ■ e Ax ^ is log- 
normally distributed, and the standard result of the log-normal distribution 
yields 

E{e Ax{ tl ) e Ax{t 2 ) . .. e Ax{t n )^ = ew h A 2 Cx (0) + A 2 Y,C x (tj - U)\. 

{ i<j J 
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□ 

Using this expression, we first have E[X(t)] = E[X(0)] = E[k e^ Xeq+x ^] = 
koexp(Px eq + ^P 2 C x (0)). For the three-step correlation function, the station- 
ary of X(t) reduces £J[AA(0)AA(*i)AA(ii + t 2 )) to 

£?[AA(0)AA(ti)AA(ti+t2)] 

= E[{X(0) - E[X(0)]}{X(h) - £7[A(ti)]}{A(ti + t 2 ) - E[X(t\ + 1 2 )}}} 
= E{X(0)X(h)X(h + 1 2 )} + 2{£[A(0)]} 3 

- E[X(0)](E{X(0)X(ti)} + E{X(0)X(h + 1 2 )} + E{X(ti)X(t\ + 1 2 )}). 

Using (A. 6), it is further simplified to 

£[AA(0)AA(ii)AA(ii + i 2 )] 

= ^e 3/3:re9+3 ' 92c ' :,:(0)/2 {e /32[c '" (tl)+c,a;( * 2)+c '" ( * 1+t2)] - e^ ^ 

_ e P 2 c x (t 2 ) _ e f3 2 c x (h+t 2 ) + 2 | 

Similarly, expanding the individual terms in the four-step correlation func- 
tion and using the stationarity and (A. 6) provide 



£[AA(0)AA(ti) AA(*i + t 2 )AX(h +t 2 + 1 3 )] 

_ u4A/3x eq +2l3 2 C x (0) 
— ft, e 

X { e / 32 [^(*l)+^( i 2)+^(*3)+C*o ; (tl+i2)+C , a; (i2+t3)+C , a ; (tl+t2+t3)] 



_ e P 2 [C x {t 1 )+C x (t 2 )+C x (t 1 +t 2 )} 
_ e /9 2 [C x {h)+C x (t 2 +t 3 )+C x (t i +t 2 +t s )] 
_ e P 2 [Cx (ti +t 2 )+C x {t 3 )+C x (h +t 2 +t 3 )] 
_ e /3 2 {C x {t 2 )+C x (t 3 )+C x (t 2 +t 3 )} , e l3 2 C x {ti) 
+ e P 2 C x (t 2 ) + e /3 2 C x (t 3 ) + e P 2 C x (tx+t 2 ) 

P 2 C x {t 2 +t 3 ) , p p 2 C x {ti+t 2 +t 3 ) 



+ e IJ'C x (t 2 +t 3 ) + e P'U x (ti+t2+t3) _ 3} 
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